
warning('off','all')
dowaveletfilter = true;

children = dir;
rootdir = cd;
piaTracking = [];
trajectoriesDDD = [];
trajectoriesD = [];
trajectoryinfo = [];
runningdata = [];
eyeDisplacement = [];
eyeRadius = [];
raw = [];
infotable = [];
objType = '25x_v2';

keyword = 'noise';

forbiddenword = 'moving';

for n = 1:length(children)
    if children(n).isdir && children(n).name(1) == 'S'        %iterate through subfolders
        cd([rootdir,'\',children(n).name])
        cursession = str2double(children(n).name(2:end));
        temp = dir('*preselected2.mat');
        temp2 = dir(['*',keyword,'*.mat']);
        temp3 = dir(['*',forbiddenword,'*.mat']);
        if isempty(temp) || isempty(temp2)
            continue;
        end
        
        %little patch to account for a mistiming issue with one directory:
        if strcmp(cd,'\\research.files.med.harvard.edu\neurobio\GU LAB\TWO-PHOTON MICROSCOPY\visualStim\thy1gcamp\013\S2')
            prestimTime = 3.5;
        else
            prestimTime = 4.5;
        end

        try
            temp = dir('*groupedResults.mat');
            if ~isempty(temp)

                regExpindicator{1} = keyword;
                regExpindicator{2} = forbiddenword;
                [a,b,c,d,e,f,~,h,i,j] = poolStimData_v5_3(regExpindicator,[],cursession,objType,prestimTime);
            else
                [a,b,c,d,e,f,~,h,i,j] = poolStimData_v5_0(keyword,[],cursession,prestimTime);     %older iteration used v1_0
            end
        catch
            disp(['unable to do session #',num2str(cursession)])
            continue;
        end
        if prestimTime == 3.5
            a = [nan(size(a,1),10),a];
            d = [nan(size(d,1),10),d];
            e = [nan(size(e,1),10),e];
            h = [nan(size(h,1),10),h];
            i = [nan(size(i,1),10),i];
        end
        infotable = [infotable;j];
        trajectoriesDDD = [trajectoriesDDD;a];
        trajectoryinfo = [trajectoryinfo;b];
        try
            piaTracking = [piaTracking;c];
        catch
            c = rmfield(c,'preStimTime');
            c = rmfield(c,'postStimTime');
            c = rmfield(c,'stimparams');
            piaTracking = [piaTracking;c];
        end
        trajectoriesD = [trajectoriesD;d];
        runningdata = [runningdata;e];
        eyeDisplacement = [eyeDisplacement;h];
        eyeRadius = [eyeRadius;i];
        try
            raw = [raw;f];
        catch
            f(:).running = [];
        end
    end
end
cd(rootdir);

if dowaveletfilter
    for n = 1:size(trajectoriesD,1)
        cur = trajectoriesD(n,:);
        missing = isnan(cur);
        if sum(missing) > 0 && sum(missing)<(.25*length(cur))
            cur = resample(cur,1:length(cur));
        elseif sum(missing)>(.25*length(cur))
            continue;
        end
        trajectoriesD(n,:) = wdenoise_NVC_v2(cur);
    end
    
end

%% calculate baseline for each vessel on a per-loop basis to reduce noise (in the baseline)

baselineframes = 25:45;      %for 4.5s baseline
[DDD,baselinemean,baselinestd,DDDz] = groupDDDcalculate(trajectoriesD,infotable,baselineframes);


%% filter out trajectories with too much fluctuation in diameter before stimulus:
bsrange = 10:45;
sloperange = 10:45;
bs = nanstd(trajectoriesD(:,1:45),[],2);        %check for standard deviation in the pre-stim period
maxslope = nan(size(bs));           %check for changes in diameter
for n = 1:size(DDD,1)
    if sum(isnan(DDD(n,sloperange))) > .45*length(sloperange)
        continue;
    end
    temp = smooth(DDD(n,sloperange),10,'rlowess');
    temp = abs(diff(temp));
    temp = temp(1:(length(temp)-5));
    maxslope(n) = max(temp);
end
bsthresh = .35;
slopethresh = .5;
ind = find(bs>bsthresh & maxslope > slopethresh);
DDDfilt = DDD; DDDfilt(ind,:) = nan;

warning('on','all')

